% solver_nopump_noLA_AOVR_MVR
% used to solve differential equations without a pump 
% aortic and mitral valve are both R-models

% initial procedure: perfectly closing mitral and aortic valve
% forward resistance = finite; backward resistance is infinite
% data are stored in resistance-matrix
% first element = forward resistance; 3rd element is backward resistance
% values are switched based on sign of PLV-PLA and PAO-PLV 
% if sign = -1: forward flow
% if sign = 1: backward flow
% resistance = resistance(sign(**)+2)
% for forward flow, the 1st element is used; for backward flow the 3rd element is used


Rmv(1)=Rmvfwd;
Rmv(2)=Rmvbwd;
Rmv(3)=Rmvbwd;
Rao(1)=Raofwd;
Rao(2)=Raobwd;
Rao(3)=Raobwd;

clear time_plot;
clear Pao_plot;
clear PLV_plot;
clear VLV_plot;
clear Qao_plot;
clear Qmv_plot;

niter=1;
diffBP=10;
SBPold=100;
DBPold=70;
DBPr2old=68;

while diffBP > prec_BP

   options=[];
   
   % There are 6 unknowns to be solved
   
   PLV=zeros(npoints,1);			% left ventricular pressure
   VLV=zeros(npoints,1);			% left ventricular volume
   Qmv=zeros(npoints,1);			% mitral flow
   Qao=zeros(npoints,1);			% aortic flow
	Pao=zeros(npoints,1);			% aortic pressure; array initialized at value DBPold
 	Pr2=zeros(npoints,1);			% presssure in 4-element WK model
    
    
   if niter==1 						% 1st iteration: assume initial values
      VLV(1)=VdeadV+Pvenous/EminV;
      Pao(1)=DBPold;
      Pr2(1)=DBPr2old;
   else   								% higher iterations: use last value of previous cycle
      HandleCheckboxFixedLVEDV=findobj(gcf,'Tag','CheckboxFixedLVEDV');
      if get(HandleCheckboxFixedLVEDV,'Value')==1;
  	      VLV(1)=VdeadV+Pvenous/EminV;
      else   
	      VLV(1)=VLV_last;
      end
      Pao(1)=Pao_last;
      Pr2(1)=Pr2_last;
   end      
      
   % the overall iteration loop should start here
   % determination of end of isovolumic contraction period
   % cross-over of isovolumic LV pressure and aortic decaying pressure
   % find the point in the interval 0-t(Emax)
   
	options = odeset('RelTol',0.0001);   
      
   t_solve=(1:npoints)*t_interval';
   [t,y]=ode45('solve_ODE',t_solve,[VLV(1) Pr2(1) Pao(1)],...
        options,R,C,L,Z0,time,E_t_LV,dEdt_LV,Rao,Rmv,VdeadV,Pvenous);
   VLV=y(:,1);
   Pr2=y(:,2);
   Pao=y(:,3);
   
   PLV=E_t_LV.*(VLV-VdeadV);
   Qao=(PLV-Pao)./Rao(sign(Pao-PLV)+2)';
	Qmv=(Pvenous-PLV)./Rmv(sign(PLV-Pvenous)+2)';

   % plot cycle results in figure 2
   
   PLV_plot((niter-1)*npoints+1:niter*npoints)=PLV;
   VLV_plot((niter-1)*npoints+1:niter*npoints)=VLV;
   Pao_plot((niter-1)*npoints+1:niter*npoints)=Pao;
   Qao_plot((niter-1)*npoints+1:niter*npoints)=Qao;
   Qmv_plot((niter-1)*npoints+1:niter*npoints)=Qmv;

   time_plot=((niter-1)*npoints+1:niter*npoints)*t_interval;
   
   axes(HandleAxesPao);
	ylabel('Pressures');
   plot(time_plot,Pao,'y',time_plot,PLV,'c');
   set(HandleAxesPao,'XLim',[0 max(time_plot)]);
   set(HandleAxesPao,'YLim',[min([min(Pao_plot) min(PLV_plot)]) max([max(Pao_plot) max(PLV_plot)])]);
   hold on;
   
   axes(HandleAxesQao);
   ylabel('aortic and mitral flow');
   plot(time_plot,Qao,'y',time_plot,Qmv,'c');
   set(HandleAxesQao,'XLim',[0 max(time_plot)]);
   set(HandleAxesQao,'YLim',[min([min(Qao_plot) min(Qmv_plot)]) max([max(Qao_plot) max(Qmv_plot)])]);
   hold on;
   
   axes(HandleAxesPLV);
   ylabel('LA and LV volume');
   plot(time_plot,VLV,'c');
   set(HandleAxesPLV,'XLim',[0 max(time_plot)]);
   set(HandleAxesPLV,'YLim',[min(VLV_plot) max(VLV_plot)]);
   
   hold on;
   pause(.01);
   
	% calculate converge criterium   
   
   if niter>2
      diffBP=abs(max(Pao)-SBPold);
   else
      diffBP=10;
   end
   SBPold=max(Pao);
   DBPold=min(Pao);
   
   VLV_last=VLV(npoints);
	Pao_last=Pao(npoints);
   Pr2_last=Pr2(npoints);
   
   niter=niter+1;
end;
